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Abstract: The miscible displacement of one incompressible fluid by another in a porous media is 
governed by a nonlinear coupled system of two partial differential equations, one of elliptic 
form for the pressure and the other of parabolic form for the concentration of the fluids. 
In this paper, the concentration is approximated by a modified method of characteristics 
(MMOC) combined with dynamic finite element spaces, while the pressure and Darcy 
velocity of the mixture are approximated simultaneously by a mixed finite element method. 
By adopting the negative-norm estimate, convergence analysis and error estimates are 
established. 
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1 Introduction 


Miscible displacement is an enhanced oil-recovery process that has attracted considerable 
attention in the petroleum industry over last 60 years. It involves injection of a solvent at cer- 
tain wells in a petroleum reservoir, with the intention of displacing resident oil to other wells for 
production. This oil may have been left behind after primary production by reservoir pressure 
or secondary production by waterflooding. The economics of the process can be precarious, 
because the chemicals it requires are expensive and the performance of the displacement is by 
no means guaranteed. Mathematically, the process, i.e., miscible displacement of one incom- 
pressible fluid by another in a porous media 2 C R? over time interval J = (0,7'] is modeled 
by a nonlinear coupled system of two partial differential equations. 

In this paper, a dynamic finite element method combined with the MMOC is applied to the 
concentration equation, while the pressure equation is approximated by a mixed finite element 
method. The main idea of dynamic finite element method is that a different number of finite 
element spaces is adopted at different time level, and the approximate solution at the current 
time is projected in the L?-norm onto the next time finite element space and make it as an 
initial value. For the theoretical analysis and practical computation of dynamic finite element 
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methods, we refer the reader to [1-5] and the references cited therein. By adopting the negative- 
norm estimate, we prove that when Mh?/At is bounded, the L?-norm is optimal, and under 
the condition that Mh*/At is bounded, optimal error estimate in energy norm is showed. 


2  Characteristic-mixed methods with dynamic finite element spaces 


The mathematical model of miscible displacement problem in petroleum reservoir simulation 
is of the following form 
(a) —V-(a(c)(Vp~7(c)Vd)) = V -u = q(z,t) TEQ, ted, 
(b) 8 +u- Ve-V.(D(z,u)Ve)=(č- c), TEN, tEJ, 
(c) u-n=D(u)Ve-n=0, rEdan, ted, 
(d) c(x,0) =co(z), TEQ, t=0. 
For detailed meanings of the parameters and some regular hypotheses of the solutions and 
coefficients in problem (1), we refer to [6-8], for example. 
Basically, the methods of characteristics are to think of the hyperbolic part of Eq (1b), 


namely ¢0c/dt+ u- Vc, as a directional derivative. Accordingly, let s denote the unit vector in 
the direction of (u1, u2, $) in Q x J, and set 


(a) = (lula)? +g)? = (ui (@)? + wala)? + o(@)?)'””. 
Then Eq (16) can be written in the following form 
Oc Bee ches 

+5, —V-(DVc) + Gc = Ge. (2) 
Let V = H(div;2)N {v -n =0 on 8N}, W = {q € L?(N) : fo ade = 0}. For a, BEV, pE W, 
and 0 € L™©(Q) define the following bilinear forms 

A(0;a,8) =(—a,8), Bla,d)=~(V-a,4). 
’ ? a(0) 3 ’ ’ F 


Then we obtain the following equivalent weak formula of problem (1) 


(a) (W3E,2) + (D(u)Ve, Vz) + (de, z) = (Gz), z € H9), 
(b) A(c;u,v) + Biv, p) = (y(¢)Vd, v), veV, (3) 
(c) Blu, $) T —(9, $), . $ € W, 


where 0 < t < T and c(z,0) = co(z). 

Suppose the time step At = T/N with tn = nAt, n = 0,1,---,N, where N is a positive 
integer. Let Va x Wn be a Raviart-Thomas space of index kn associating with a quasi-regular 
triangulation of Q, such that the elements have diameters bounded by hy in every time tn. 
Define V, = {v € Vnjv:n = 0 on NQ}, Wn = W,/{¢ = constant on 2}. It is clear that 
Va x Wn C V x W. Next let Mn C WŁ (NQ) be a standard finite element space for a Galerkin 
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method of index J, associated with another quasi-regular polygonalization of Q, such that the 
elements have diameters bounded by h?. Set 


k = info<n<nkn, l=infocngnln, he =SUPpcn<n hes Rp = 8UPo<n<n hp: 


The reader should not be confused with h?(h?) and hf(hk), where the former is the spatial 
step size at time tn, and the latter means h,(hp) to the kth power. 

At each time level n, we assume {c”, u”, p”} is approximated by {C",U", P™} € Mn x Vn x 
Wn, respectively, and for the initial approximation, we assume C° is the elliptic projection of 
co(x) onto Mo. Now we define our approximation scheme as follows. 

1) Find Ĉĉ” € Mmn+1, such that 


(@C™, z) = (@C”, 2), z E€ Mnii, n = 0,1,- „N-41. (4) 


2) When Ĉ” is known, we get {0”+t1, U” +}, Pr+1} € Mn+1 X Vn+1 X Wn+1, such that for 
z € Mn41,V E Vagi and w € Wn41 


corti — on n+1 +1 n+l n+l getignti, 
(Sa) + (DUVE, Ve) + (Ho, 2) = G z) (5) 
A(C";U"*}, v) + B(v, P?*1) = (y(")Va, v), (6) 
B(U"™*!,w) = -(q""1,w), n=0,1,---,N-1, (7) 


where yee 
ae) At). 
$(z) 
Noting that from (4), we can easily see that when different finite element mesh or interpo- 
lation functions are used at time level t = tn and t = tn41, we must project the approximate 


Č” = 6™(4) = Ô” (z - 


solution C” to the next time finite element space Mn+1, and make it as an initial value. In 
practice, once we know C” by (4), we can obtain {U"+!, P"+1} by (6)-(7), and then get C?*+4 


by (5). 
For convenience of theoretical analysis in the following section, we introduce the following 
two useful projections 


A(c(t); Rault), v) + Biv, Rap(t)) = (y(e(t)) Vd, v), vEVn, ted, 
B(Rnu(t), w) = —(q(t),w), WE Wn, te J. 


(8) 


(D(u(t))V(Rne(t)), Vz) + (Rne(t), z) + (Q(t) Rne(t); z) 
= (D(u(t)) Velt), Vz) + (c(t), z) + (@(e(E), z) 


= (E0), 2) = (ult) + Velt),2) + (clt),2) + (GA), 2), 7EMn, tE. (9) 


For the projection solutions defined in (8)-(9), we have the following approximate results. 
Lemma 1"! There exist a constant K independent of hp and he, such that 


lu — Rrul|r(s.r(aivy) + lip- Raple) < K (lel noc a;e+1(aivy) + lpllregm ABT, 


le — Rell z-( 3:12) + helle — Rnelln=(y42) < Kllelrægemyhe 
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3 Main results 


In this section, we shall show that the optimal-order convergence rate rests with the total 
times of changing finite element spaces and discrete parameters. Set 


e” =U" — Rau", 0° = C” — Rac”, p° = u” — Ryu”, AX = c — Rac”, n=0,1, ,N, 
A” = c” — Rayic”, 6" = C” — Rapi”, n=0,1,---,N—-1. 


Lemma 2 Suppose that {C,U, P} are the solutions of (4)-(7) and {Rc, Ru, Rp} are defined 
by (8)-(9), respectively. Then 


jurt _ Ropu” t! |\y y prt: — RntiPniillw < Kjatt pa Ĉĉ”). (10) 


Lemma 3 Suppose that {C,U, P} are the solutions of (4)-(7) and {Rc, Ru, Rp} are defined 
by (8)-(9), respectively. Then 


|p2an+2 || a \|¢20" ||? + At(D(U™*1) Vet, vont) 
< K{At(jor+" ||? + JÊ" ||?) + En}, (11) 
where 


= At? ( 


els 


|S sal L?( (tnstn+1; L?) +3 ee Aas) 


L2 (tnitn+i; H!) 
HAHAY + Jelic gr + At RH elc mn) 
+At h+? (lul|Z (J; e+ (div) + lpllZ (s;m) 


Proof It follows from (5) and (9) as well as the L?-projection (4) that 


grti _—6 wen ane 
(¢——.z) + (DUH), v2) 
= dc nri n+l n+l bial et n+l antign+i 
= (jo a + Ut. Ve |- 2) - (APH: z) — (@rttontt, z) 


+([u"t? — U+]. Yet, z) + ([D(u"+!) — D(U”+)] - V Rasic” tt, Vz) 


e = x2) Š (oe 


yn n \n n gn 
At = 2) 3 (0 a 2) 7 (° a 12); (12) 

Let z = 0"*1 be the test function, for the estimates of the right-hand side of (12), we should 
only pay special attention to the estimate of the seventh term. Since if M, = Mn41, i.e., the 
same finite element spaces between the time levels tn and tn41 are used, then it follows from (4) 
that à” = Âr, and so that this term reduces to zero. Otherwise, by utilizing the negative-norm 
estimate, we have 


(oro) 


lA 


_ lei 


IA 


K(A hE Hell Eec) + elo Ni. (13) 
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Then we conclude by Lemmas 1,2, and the well known estimates for MMOC that 


aa [letet — eter iP] + (DU) ve"4, ver) 


< K\jor*4||? + Ki 6" |? + KAt( 


ž Zl ) 
OT? IlL? ta,tny1;2?) L2 (tnitn41;L?) 


+K (At) nz] 2 


t L2 (tr jtng13;H'+1) 


+ K(At)-7h2"*4 lelle) 


+K Ylle egmay + KAZ? (lull? ecs acai + lole aa) 
+e] Vor}. (14) 
Multiplying (14) by 2At and choosing e€ sufficiently small, we then complete the proof. 


Lemma 4 Let M be the total times of changing finite element spaces for concentration in 
the time direction, then we have 


N-1 
JON? + So Advert? 
n=0 
< 2 | 2 214-2 
mt ( Or? İİ LZ(J;L2) +| Ot l'n2(7; gaa) t Hele Is L2(J;Hî+1) 
+K ABE? (Jul (ggr (div)) + lllz aa) 
Mh? 
+K (= Paridis, HI+1(Q9))° (15) 


Proof It follows immediately from Lemma 3 and the projective relation (4). 

Combing the estimate in Lemma 4 with that in Lemma 1, it is easy to get the following 
main result. 

Theorem 1 Suppose {c,u,p} and {C,U, P} are the solutions of (1) and (4)-(7), respec- 
tively. Under some regularity assumptions on the solutions and coefficients, we have the fol- 
lowing error estimates for the concentration that 


_—ony2 — 2 Mh? 21+2 2k+2 

pe OC ERO i ag E (16) 
3 At||V(c" — o)? = o(a? + (bs + 1) ne! + he) (17) 
= At coe EBE P 


The size of At term depends principally on Is] and Iê. The spatial terms depends prin- 
cipally on the H'+! and H+! norms of c and u, p 

Theorem 2 Under the assumptions of Theorem 1 and the extra assumption that c € 
W21 (J; L?), the errors in velocity and pressure are bounded by 


m Mh? 
max, (llu” — Ul + lp" — P”) = o(a $ (= + ijae + pane: (18) 
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